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Abstract 

Exact numerical methods and stochastic simulation methods are developed to study the force 
stretching single RNA issue on the secondary structure level in equilibrium. By computing the 
force-extension curves on the constant force and the constant extension ensembles, we find the two 
independent methods agree with each other quite well. To show the precision of our methods in 
predicting unfolding experiments, the unfolding forces of different RNA molecules under different 
experimental conditions are calculated. We find that the ionic corrections on the RNA free energies 
alone might not account for the apparent differences between the theoretical calculations and the 
experimental data; an ionic correction to the persistent length of single-stranded RNA should be 
necessary. 
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I. INTRODUCTION 

In the past few years, enormous theoretical efforts have been devoted to understanding 
folding/unfolding phenomena of proteins and RNA observed in single-molecule experiments. 
Diverse methods including molecular dynamics □ D , Monte Carlo method □ □ □ , and 
other theoretical models 6[ have been developed. However, these studies mainly focused 
on the dynamical behavior of proteins under force, and few concerned about RNA j5|. To 
fill this gap, we recently developed kinetic Monte Carlo simulation methods to investigate 



the RNA kinetic 
structure level 



behaviors in constant force and constant extension ensembles on secondary 
1 

[. In addition to the intriguing nonequilibrium phenomena, the most 
direct application of our simulation methods is to investigate the relative simple unfolding 
behaviors in equilibrium 0,0]. Different from complex protein unfolding behaviors even in 
equilibrium state [3], force unfolding RNA has been showed to be solved in an exact numerical 
way in the constant extension ensemble [10], though the extension to the constant force 
ensemble were not reported till now. One of natural question arises whether our simulations 
are accurate enough comparing to the numerical method. In this report, we address this 
question. 

The organization of this paper is as follows. We first ,In Sec. ITT1 simply review the Monte 
Carlo methods developed by us. Then the exact numerical methods for the constant force 
and the constant extension ensembles are showed. In Sec. II I II we compare the simulation 
and numerical methods in the two ensembles. We particularly point out the importance 
of persistent length of RNA in predicting unfolding forces. Finally Sec. I1VI presents our 
conclusion. 



II. THE MODEL AND METHODS 

According to the difference of the external controlled parameters, the RNA unfolding 
experiments can be carried out under constant extension and constant force, i.e., the constant 
extension and the constant force ensembles [llj|. one of apparatus for the constant extension 
ensemble is sketched in Fig. ^ a single RNA molecule is attached between two beads with 
RNA:DNA hybrid (double-stranded DNA or dsDNA) handle; one bead is held by a pipette, 
and the other is in a laser light trap. In practice, although two identical handles connect 
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the RNA, only one handle is considered in order to simplify our theoretical calculation; It 
should not change following discussions. By moving the position of the pipette, the distance 
between the two beads and the force acting on the bead in the light trap can be measured 
with high resolution. On the contrary, constant force can be imposed on the RNA molecules 
with feedback-stabilized optical tweezers capable of maintaining a preset force by moving 
the beads closer or further apart. 

A. Monte Carlo simulation methods 

The Monte Carlo algorithm is simply reviewed in this section. 

n 

First is the method for the constant extension ensemble Two simplifications have 
been made in our model. We suppose that changes of the extensions of RNA and the handle 
proceed along one direction. Physical effects of the beads are neglected. Consequently, any 
state of the system can be specified with three independent quantities, the position of the 
bead with respect to the center of the optical trap, x tw , the end-to-end distance of the handle, 
x ds , and the RNA secondary structure S, i.e. the system in z-state (Si, xf", x ds ). Here we 
do not include x ss , the extension of the RNA for the sum of individual extensions satisfies 
constraint condition, x = x tw + x ds + x ss , where x is the distance between the centers of the 
light trap and the bead held by the pipet, and it also is the external controlled parameter 
in the constant extension ensemble. The move set for this system is as follow, 

(C „tw ds\ I o tw ds\ ■ / • 

(S i ,x?,x* s )^(S i ,xr±8,x ds TS), (1) 
(S i ,xf,x ds )^(S i ,xf,x ds ±5). 

Unfolding the single RNA for the constant extension ensemble proceeds in an extended 
conformational space C(l) x R tw x R ds , where C(l) is the RNA secondary structural folding 
space, R tw = (0, +00), R ds = (0, Z<f s ), and ld s is the contour length of the dsDNA handle. 
Given the system state i, its whole energy is 

Ei(x) = AG° + W tw (x^ ) + W ds (x ds ) + W ss (x° s , m), (2) 

where is the free energy obtained from folding the RNA sequence into the secondary 
structure Si, and the elastic energies of the optical trap, the handle, and the single-stranded 
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part of the RNA are 



W tw (xf) = l -k tw xf\ 

,.ds 

* 

rds/ds\ I e Jj 



W d \xf) = / f ds (x r )dx r , (3) 







W ss (xr,n t ) = xi'fixi^m) -I ' ' x. a (f,ni)df, 



o 



respectively. In the expression W , fds(x r ) is the average force of the handle at given 
extension x', 



k B T ( 1 la;' 



where is the persistence length, respectively. In the expression W ss , x ss (f',ni) is the 
average extension of the single stranded part of the RNA whose bases (exterior bases) is rii 
at given force /', 

x ss (f',ni) = n;6 ss [coth(|-^) - (5) 
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where b ss and l ss are the monomer distance and the Kuhn length of the single-stranded 
RNA, respectively Q, Q|- Note that f(xf,rii) is the inverse function of x ss (f',ni). 

Then is the simulation method for the constant force ensemble Q . We proposed an energy 
expression on the coarse-grain level for the given secondary structure Si under constant force 

/, 

3t(/)wAG?-n,X0(/), (6) 

where g(f) = kBTb ss /l ss \nsmh{u)/u and u = l ss f ' jk^T. In contrast to the constant ex- 
tension ensemble, the RNA secondary structure S can completely specify any state of the 
constant force ensemble. Therefore, the move set for this ensemble is the same with the set 
for RNA folding without force, i.e., its unfolding space is C(l). 

Given the move sets and the unfolding conformational spaces, the RNA unfolding for the 
two ensembles can be modelled as a Markov process in their respective spaces. Define the 
transition probabilities kij from i-state to j-state satisfying kij = r" 1 exp(— A£ , j J /2/csT), or 
the symmetric rule (l^ . where AEij = Ej — E i7 t q is used to scale time axis of the unfolding 
process. We use a continuous time Monte Carlo algorithm to simulate unfolding process 

he; 



The measurement quantities (A) for the two ensembles can be calculated by (A)(n) = 
Ai(U + i — ti), where Ai is the A-value in state i, and t{ is the inner time of the Monte 
Carlo simulations. For the constant extension ensemble, A could be the force exserted on 
the bead, / = k tw x tw in the light trap, or the bead-to-bead distance x bb = x ds + x ss . While 
for the constant force ensemble, A is the molecular extension x under the constant force /, 
and X{ = x ss (f,rii). The simulation time is 2 x 10 6 r o . 



B. the exact numerical methods 



Compared to difficult protein folding prediction, the RNA secondary structure prediction 

n 

has achieved great success (18[ . In particular the partition function method developed later 
provided strongly physical foundation Q|. Recently, this method was generalized to the 
case of RNA unfolding in the constant extension ensemble Q| • In the present work, we are 
not ready to choose the formulae presented in Ref. In addition to be consistent with the 
formulae for the Monte Carlo simulations, the complicated polymer model of single-stranded 
DNA (ssDNA) therein might not result in many advances in predicting and understanding 
the RNA unfolding phenomena. 

The key idea of the partition function method is that the partition function over all 
secondary structures of a given RNA can be calculated by dynamic programming. Given 
the partition function Q(i,j,n) on the sequence segment [i,j] with exterior bases n, its 
recursion formula is as follows, 

Q(i,j,n) = I6k,j-i+i + qb(i, A+j — n) 

j-l k-i+l 

+E E ( 7 ) 

k=i m=l 

xqb(k + 1 , m + A + j — n) , 

where the partition function qb(i,j) on the sequence segment [i,j] for which the i and j bases 
are paired; Vienna package 1.4 provides their calculation codes j^ . 

For a given RNA sequence consisting of iV nucleotides, define the total partition function 
for the constant extension and the constant force ensemble Z^{x) and Z^(f), respectively. 
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According the energies mentioned in last section, their expressions can be written as 

N fids rnb ss 

Z N (x) =J2 dxdS / dx ss Q(l,N,n) 

xexp{-(3E{x,x ds ,x ss ,n)) (8) 

and 

N 

ZnU) = E Q0-> N > n ) ex p(-^(/' n )) ( 9 ) 



where the elastic energy E(x, x ds , x ss , n) = W tw (x — x ds — x ss ) + W ds (x ds ) + W ss (x ss , n) and 
E(f, n) = nx g(f). Correspondingly, the measurement quantities for the constant extension 
ensemble are the average force (/) = —kBTdZ^(x)/dx and the average extension (x bb ) = 
x — (f)/k tw , and (x) = k B TdZ N (f)/df for the constant extension ensemble, respectively. 

III. COMPARISON OF THE EXACT AND SIMULATION METHODS 

To compare the exact and simulation methods discussed above, we calculate extension- 
force curves of three small RNA, p5ab, p5abcAA and p5abc in equilibrium. Their native 
states under experimental condition are showed in Fig. ^ These molecules have been 



studied by the experiment 



llj and simulation j?l Isj]. We first choose the widely used 



parameters for our computation: temperature T = 298K, b ss = 0.56 nm, l ss = 1.5 nm, 



P ds = 53 nm 



12L Q], Ids = 320 nm, and k tw = 0.2 pN/nm and the free energy 



parameters for RNA secondary structures at standard salt concentrations: [iVa + ] = 1M 



and [Mg 2+ ] = 0M |20|. 

Fig. El shows these extension-force curves for the sequences for the two ensembles. We 
find that the two independent methods achieve highly consistence. In particular, the three 
curves of the molecules for the constant extension ensemble also agree with the experimental 
measurements very well in quantity: the extensional transition of P5ab are all-or-none; 



while P5abc has an intermediate state Interestingly, we note that, P5abcAv4 although 
has been observed as two-state molecule in the experiment, a weaker intermediate state 
presents in the constant extension ensemble, while it cannot be observed in the constant 
force ensemble. 

If we purchase the precision of our methods, quantitative comparison between the theoret- 
ical molecular unfolding force f u and the experimental measurements of course is essential. 
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But we find that they do not coincide: in the constant extension ensemble, the experimental 
unfolding force of P5ab is 13.3 pN, of p5abcAA is 11.4 pN. and of P5abc is 8 pN [llj]; while 
our calculations are 18.4 pN, 15.8 pN and 12.2 pN, respectively. So what causes result to 
the larger differences between the experiment and the theory? The experiment and previous 
theoretical works contributed the differences to the change of free energy of the RNA sec- 
ondary structure; this change results from the different ionic concentration of the experiment 
and the standard condition: in the RNA unfolding experiment, Na + = 250 mM and with 



and without Mg 2+ = 10 mM 11, 21]. To reproduce the experimental ionic condition, 
a correction on the energy of a base pair equal to — 0.193fcsT ln([iVa + ] + S^Mg 2 " 1- ] 1 / 2 ) has 

ecules from 



been applied j2l|. Their values are summarized in Tab. HI Besides the three mo 
Ref. other unfolding forces of the molecules published in the lectures 22, 



23j are also 



listed there. We still see that the ionic correction cannot explain the derivation between the 
theory and experiment. 

Considering that the free energy parameters of the RNA secondary structure were mea- 
sured in bulk experiments, one might doubt whether they can be used in single-molecule 
studies as well as we thought before jsj]. On the other hand, however, it is known that the 
mechanical parameter, the persistent length l ss is also sensitive to ionic condition. Although 
this parameter indeed were measured under a similar experimental conditions with the small 
RNA unfolding experiment (see Ref. 12 1), their validity for describing small molecules is 
questionable. Recent FRET experiment measured that for shorter ssDNA l ss is about 2.2 
nm at Na + = 250 mM 24J. If we choose this value in our calculation, the predicted unfold- 
ing forces are closer with the experimental measurements; see Tab. d Of course, we cannot 
exclude the intrinsic limitation of our coarse-grain model. For example, another possible 
force work formula has been used in the constant force ensemble j?|. 



IV. CONCLUSION 



In this work, we review the Monte Carlo methods and develop the exact numerical meth- 
ods to study the force stretching single RNA molecules issue. We respectively compare the 
two independent method in the constant force and extension ensembles, and find that they 
agree with each other quite well. We also point out that only ionic correction on the RNA 
secondary structure alone cannot explain the larger discrepancies of the unfolding forces 
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between the theoretical prediction and the experimental measurement; the ionic correction 
on the RNA molecular mechanical properties should be important. 

Although the results of the exact numerical method are consistent with the Monte Carlo 
method when force stretches single RNA in equilibrium, it does not mean the former can 
completely replace the later. Such situation is similar with the study of 2-dimension Ising 
model in condense matter physics Compared to the exact method, the Monte Carlo 
method would be more sophisticate in dealing possible more complicated experimental con- 
dition. For instance, recent simulation work could include pseduknots structure j^l, while 
the exact partition function technique would be hardly to realize. In our point of view, 
Monte Carlo simulation is more important in studying single molecular non-equilibrium be- 
havior produced by mechanical force, such as folding/ unfolding trajectories, force- hysteresis 
phenomena and unfolding force dependance on loading rates etc p, Q . 

We thank Professor H.-W. Peng for many helpful discussions in this work. 
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TABLE I: The unfolding forces f u of different molecules under different experimental conditions. 
The experimental data are from the previously published data 
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22 



231 ] . The theoretical values 



are from the exact numerical methods developed above, where /* , i = 1,2,3 represent the unfolding 
forces without the ionic correction, with the ionic correction on the free energy and with the ionic 
and the persistent length corrections, respectively. Here We do not show the P5abc unfolding force 
for it is not reversible in Mg 2+ due to the presence of tertiary interactions. 



Molecule 


temperature (K) Na 


+ (mM) Mg 2 


+ ( 


mM) fl (pN) / 


I (PN) / 


, 3 (PN) f u xp (pN) 


P5abc 


298 


250 





12.2 


11.4 


10.0 


7.0-11.0 


poly(dA-dU) 


293 


150 





12.3 


11.0 


9.3 


9.0 


P5abcAA 


298 


250 





15.8 


14.8 


13.2 


11.4 ±0.5 


P5abcAA 


298 


250 


10 




15.4 


13.8 


12.7 ±0.3 


P5ab 


298 


250 





18.4 


17.4 


15.7 


13.3 ± 1.0 


P5ab 


298 


250 


10 




18.0 


16.2 


14.5 ± 1.0 


CG hairpin 


293 


150 





25.8 


24.4 


22.4 


17.0 


poly(dC-dG) 


293 


150 





25.1 


23.8 


21.7 


20.0 
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FIG. 1: Theoretical model and RNA sequences and their native structures studied in present work. 
The structures are folded by Vienna RNA package 1.4. The equilibrium and kinetic behaviors of 
these three RNAs, p5ab, p5abcAA, and p5abc have been studied in detail jll| . 
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FIG. 2: Comparison of the exact and simulation force-extension curves in equilibrium for P5ab, 
P5abcAA and P5abc in the two ensembles. The different symbols are from the simulation methods, 
and the different lines are from the exact methods. They agree with each other very well. 
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